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ABSTRACT 

The Padoan and Nordlund model of the stellar initial mass function (IMF) is derived from low order 
statistics of supersonic turbulence, neglecting gravity (e.g. gravitational fragmentation, accretion and 
merging) . In this work the predictions of that model are tested using the largest numerical experiments 
of supersonic hydrodynamic (HD) and magneto-hydrodynamic (MHD) turbulence to date (~ 1000 3 
computational zones) and three different codes (Enzo, Zeus and the Stagger Code). The model 
predicts a power law distribution for large masses, related to the turbulence energy power spectrum 
slope, and the shock jump conditions. This power law mass distribution is confirmed by the numerical 
experiments. The model also predicts a sharp difference between the HD and MHD regimes, which is 
recovered in the experiments as well, implying that the magnetic field, even below energy equipartition 
on the large scale, is a crucial component of the process of turbulent fragmentation. These results 
suggest that the stellar IMF of primordial stars may differ from that in later epochs of star formation, 
due to differences in both gas temperature and magnetic field strength. In particular, we find that 
the IMF of primordial stars born in turbulent clouds may be narrowly peaked around a mass of order 
10 M , as long as the column density of such clouds is not much in excess of 10 22 cm~ 2 . 
Subject headings: ISM: kinematics and dynamics — stars: formation — turbulence 



1. INTRODUCTION 

In the turbulent fragmentation model of Padoan & 
Nordlund (2002), the mass distribution of gravitationally 
unstable cores in turbulent clouds is derived from low 
order statistics of supersonic turbulence (the one-point 
pdf of gas density, and the first-order scaling of velocity 
differences -though we refer to the velocity power spec- 
trum, as a second order proxy of the first order scaling) 
and from focusing on the fundamental flow geometry of 
shock compressions. This is a natural (if not traditional) 
approach to the solution of a very complex turbulence 
problem, because a simple solution must be based on 
low order statistics giving the basic variance and scaling 
of the process. But low order statistics cannot capture 
the flow geometry. Instead of relying on closure mod- 
els of a statistical nature, an alternative is to impose 
some knowledge of a fundamental flow feature that would 
otherwise be hard to capture with statistical quantities. 
This is especially true because of the strong intermittent 
nature of turbulent flows, meaning that the most im- 
portant flow structures could be properly accounted for 
only by very high order statistics, so relatively low order 
statistical closure models are bound to largely overlook 
those important structures. A phenomenological model 
centered on the geometry of those structures, in combi- 
nation with low order statistics is therefore a valid alter- 
native, and it is the approach of choice in the work of 
Padoan & Nordlund (2002). 

The model assumes that: i) the turbulence has a power 
law velocity power spectrum (the model really uses the 
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first order scaling of velocity differences, but that is re- 
lated to the power spectrum anyway); ii) cores are formed 
by shocks in the turbulent flow and have size and density 
that scale as the postshock layer thickness and density; 
iii) the number of such shocks scales self-similarly as the 
inverse of the cube of their size; iv) the condition for the 
collapse of small cores is that their mass exceeds their 
Bonnor-Ebert mass, derived from the lognormal proba- 
bility density function (pdf) of the gas density indepen- 
dently of the core mass. 

These assumptions are all reasonable for an approx- 
imately isothermal and supersonic turbulent gas. The 
first assumption is a well established result for turbulent 
flows. The second assumption is suggested by the simu- 
lations, where the unstable cores are always found to be 
the densest regions of postshock sheets or filaments. The 
third assumption is motivated by the very large Reynolds 
number of the turbulence in star-forming clouds, which 
is expected to generate a very extended inertial range 
of scales and possibly a self-similar flow responsible for 
the network of shocks (this self-similarity does not imply 
necessarily a hierarchical, space-filling nested structure, 
where smaller compressing regions would always be in- 
side larger ones). Finally, the fourth assumption stems 
from the fact that most of the densest gas is located 
within dense cores, so the high density tail of the gas 
density pdf should be well represented by the distribu- 
tion of the core mean density. 

These assumptions result in a power law mass distribu- 
tion reflecting the scale-free nature of the turbulence, and 
a turnover at small masses, where gas pressure competes 
with self-gravity. After integrating over the probability 
of exceeding the Bonnor-Ebert mass, the mass distribu- 
tion is given by the following formula: 
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Fig. 1. — Logarithm of projected density from a snapshot of the 
Stagger-Code HD run. 



where the mass m is in units of the average Bonnor-Ebert 
mass, 

™b E .o = 3. 3 M s (^)-' /2 (^) 3/2 , (2) 

a is the standard deviation of the gas density pdf (as- 
sumed to be a lognormal) related to the rms Mach num- 
ber of the turbulence (the sonic or the Alfvenic Mach 
number in the HD or MHD regime respectively, see be- 
low): 

a = v/ln(l + Af 2 /4) (3) 

(from here on the subscript o denotes quantities averaged 
on the outer scale of the turbulence). The coefficient C 
depends on the physical parameters and is not discussed 
here because the normalization of the mass distribution 
is not required for this work. The power law slope, x, is 
determined by the power law slope of the energy spec- 
trum, (3 {(3 w 5/3 in incompressible turbulence and (3 — 2 
in Burgers zero-pressure model), and by the shock jump 
conditions: 

x = 3/(4-/3) (4) 

for B > B CI (MHD jump conditions), and 

1 = 3/(5-2/3) (5) 

for B < B CI (isothermal HD jump conditions). The crit- 
ical magnetic field value that separates the two regimes 
is given by the condition that the postshock gas pressure 
is of order the postshock magnetic pressure, correspond- 
ing to an rms Alfvenic Mach number, Ma, of the or- 
der of the ratio of the mean gas and magnetic pressures, 
Ma ~ P g /P m . This condition gives 

Because the Galactic magnetic field strength is locally 




Fig. 2. — Logarithm of projected density from a snapshot of the 
Stagger-Code MHD run. 



6±2 fiG (Beck 2001; Han, Ferriere, & Manchester 2004), 
and perhaps larger in molecular cloud cores (Crutcher 
1999; Bourke et al. 2001), current star formation in the 
galactic disk occurs under MHD conditions. For a value 
of (3 = 1.9, found from the least dissipative simulations 
in this work, the predicted slope of the mass distribu- 
tion of prestellar cores is then x — 1.4, practically the 
same as the Salpeter slope of the stellar IMF (x = 1.35, 
Salpeter 1955). For very weak magnetic fields, perhaps 
in protogalaxies at very large redshifts, the slope would 
be x = 2.5, assuming again (3 = 1.9. Furthermore, con- 
ditions at high redshifts and very low or zero mctallicity 
would also include a larger temperature, T > 100 K (e.g. 
Palla, Salpeter, & Stahler 1983; Abel, Bryan & Norman 
2000) . The larger temperature results in a value of B CT at 
least 10 times larger than in present day environments 
with the same rms velocity and mean density, making 
the HD regime, and hence the steeper IMF, even more 
likely to occur for stellar populations at high redshift. 
At the same time, the larger temperature also shifts the 
peak of the IMF to larger masses, according to equations 
((T|) and ((2|). Although the first population III stars are 
usually believed to form in isolation in the cores of very 
small halos, population III stars formed somewhat later 
and in somewhat more massive halos (Jimenez & Haiman 
2006; Iliev et al. 2006) may indeed emerge from turbu- 
lent star-forming environments more akin to current star 
formation sites. Thus, the steeper IMF for massive stars 
corresponding to the HD regime may also apply to bona- 
fide population III stars, as well as to early population 
II stars. 

The peak of the distribution shifts to smaller masses 
with increasing Mach number and gas density, which also 
increases the abundance of low mass stars and brown 
dwarfs. For reasonable values of the physical parame- 
ters, this mass distribution becomes essentially the same 
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Fig. 3. — Compensated power spectra of the main four runs and 
least square fits in the range of wavenumbers 3 < k < 20. 



as the observed stellar IMF in Chabrier (2003). This 
coincidence suggests that the process of turbulent frag- 
mentation may play a major role in the origin of the 
stellar IMF (Padoan & Nordlund 2002), with only minor 
effects due to gravitational fragmentation, accretion or 
merging. 

A numerical estimate of the dependence of the peak 
mass on the rms Mach number gives: 



ITlr 



3.0M n - 1 - 1 M B E,o 



assuming x = 1.4 (MHD regime), and 



ITlr 



3.9M^ 1 - 7 Mbb,o 



(7) 



(8) 



assuming x — 2.5 (HD regime). In the following we relate 
m p to the peak of the stellar IMF. Since the latter can 
be accurately determined observationally only for stel- 
lar clusters, presumably formed in gravitationally bound 
cores, we may assume the gas velocity dispersion is of or- 
der the virial velocity. With this assumption, and using 
the rms Alfvenic Mach number, equation ([7]) gives: 

m p = O^MoA^XiW^io ( 9 ) 

where -/V 0) 22j -So.icb n o,4j an d T QtW are the mean column 
density, magnetic field, particle density and temperature 
in units of 10 22 cm" 2 , 10 /xG, 10 4 cm -3 , and 10 K re- 
spectively, which are typical values for star-forming cloud 
cores. This value of m p provides an estimate of the peak 
of the stellar IMF (although it should be somewhat larger 
than the IMF peak because a fraction of the core mass is 
not accreted onto the star and because some cores may 
result in binaries or multiple systems). 

The observed IMF (of multiple systems) peaks at 
roughly 0.2 M (Chabrier 2003), with no clear evidence 
of a strong departure from this value in any star-forming 
region, with the exception of the Taurus molecular cloud 
complex (Luhman 2004). As both temperature and col- 
umn density don't have very large variations from cloud 
to cloud, our fragmentation model would predict a nearly 
constant value of m p if Bo oc «g' 45 , for which there 
is some observational and theoretical support (Myers 
& Goodman 1988; Crutcher 1999; Padoan k, Nordlund 
1999; Basu 2000). Sites of massive star formation may 
have larger mean column density, but they also have 



Fig. 4. — Mass distributions of gravitationally unstable cores 
above 1 Mq, for the main four experiments scaled to a mean den- 
sity of 10 4 cm -3 , a box size of 6 pc, and a clumpfind density 
resolution / = 8%. The dashed lines show the power law derived 
from the power spectrum slope and the shock jump conditions of 
the corresponding simulations, according to the turbulent fragmen- 
tation model. The histograms are arbitrarily offset in the vertical 
direction for clarity. 



larger temperature than regions of low-mass star forma- 
tion (e.g. Ossenkopf, Trojan, & Stutzki 2001; Beuther et 
al. 2006), which may result in a similar value of m p . 

In the HD regime, assuming again virial velocity dis- 
persion and using the sonic rms Mach number, the peak 
mass given by equation ((SJ) is: 



to p = O.16M iV ( 



1.7 0.36 



^2.35 



0.22 



0,4 J 0,10 



(10) 



Assuming for example that the earliest turbulent 
conditions arc found within halos of approximately 
10 s Mq, cooled to the background temperature by non- 
equilibrium formation of HD (Johnson & Bromm 2006), 
approximate values of mean density, temperature, and 
column density are n$ = \ cm" 3 , To = 43 K, and 
N = 10 21 cm" 2 , giving m p = 10 M Q . It has been of- 
ten proposed that the IMF of population III stars should 
contain an excess of massive stars relative to the present 
day IMF. This may indeed occur with our IMF in the HD 
regime, despite its steep slope, due the high gas temper- 
ature of primordial gas, shifting the peak of the IMF to- 
ward large masses. Because of its steep slope, though, the 
IMF of primordial stars would then be narrowly peaked 
around a large mass of order 10 Mq. On the other hand, 
a larger column density of order 10 22 cm" 2 would give 
m p = 0.2 M Q , similar to the IMF peak in regions of 
present-day star formation. This example shows that 
the typical mass of primordial stars born in turbulent 
clouds cannot be estimated without a reliable value of 
the average column density of such clouds. 

2. THE SIMULATIONS 

As the turbulent fragmentation model outlined in the 
previous section neglects gravity (apart from the selec- 
tion of unstable cores), its predictions can be tested with 
numerical simulations of supersonic MHD turbulence. 
However, because the model relies on the scale-free na- 
ture of the turbulence, it cannot be easily tested unless an 
extended inertial range of the turbulence is generated in 
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Fig. 5. — Mass distributions of gravitationally unstable cores 
from the HD and MHD regimes (Enzo and Stagger Code respec- 
tively), computed with / = 16% and assuming a mean gas density 
of 10 4 cm -3 . Each mass distribution is the result of matching two 
mass distributions, computed for a computational box size of 1 pc 
and 6 pc. The Chabrier IMF (Chabrier 2003) and the fragmenta- 
tion model predictions for the power law mass distributions of each 
run are also plotted. 



the simulations, which requires both large numerical res- 
olution (large computers) and low numerical diffusivity 
(good codes). The main comparison between the MHD 
and HD regimes is here based on simulations with the 
Stagger Code, on a numerical mesh of 1, 000 3 computa- 
tional zones. In order to rule out numerical artifacts of a 
specific code, we also simulated the HD regime with the 
Enzo code (Norman & Bryan 1999), and the MHD regime 
with the Zeus code (Stone & Norman 1992), in both cases 
on a numerical mesh of 1, 024 3 computational zones. The 
three codes are quite different from each other. Enzo is 
based on a Rieman solver and PPM reconstruction, the 
Stagger code is based on a high order (5th order in space, 
4th order in time) finite different scheme, and Zeus is a 
low order (2nd order in space, 1st order in time) finite 
difference code. 

All simulations make use of periodic boundary condi- 
tions, isothermal equation of state, random forcing in 
Fourier space at wavenumbers l<fc<2(fc = l cor- 
respond to the computational box size), uniform initial 
density and magnetic field (in the MHD runs), random 
initial velocity field with power only at wavenumbers 
1 < k < 2. The simulations arc run for several dy- 
namical times to relax the turbulence at rms Mach num- 
bers of 6 or 10, before being analyzed. All results in 
this paper are for rms Mach number 10, unless other- 
wise specified . In Figures 1 and 2 we show two projec- 
tions of the density field, from the HD and MHD Stagger 
Code runs. The density field in the HD run appears to 
be significantly more fragmented than its MHD coun- 
terpart. There are two main reasons for this difference: 
i) The density contrast in the HD shocks is larger than 
in the MHD shocks, creating thinner postshock layers 
from shocks with equal sonic Mach number; ii) the HD 
postshock layers are Kelvin-Helmholtz unstable, due to 
the strong shear flow that originates in oblique shocks, 
while in most of the MHD layers the same instability 
is suppressed by the magnetic field that is amplified in 



1000 






■ i i | 




i 










— Lj^s 








100 
























mhd ; 
\ s Zeus . 


10 




Chobrier 

m 




HD \l 


1 \ 




1 




. . . 1 


. . . i 


PPM I 


. . A i 






0.1 


1.0 




10.0 





m [M ] 

Fig. 6. — Same as Figure 5, but for the Zeus MHD run instead 
of the Stagger-Code MHD run. Notice how the steeper Zeus mass 
distribution is recovered. The model predicts the Zeus mass distri- 
bution to be steeper than the Stagger-Code one, as a result of the 
steeper turbulence power spectrum in the more diffusive Zeus run. 



the compression. The turbulent fragmentation model of 
Padoan and Nordlund (2002) makes explicit use of the 
shock jump conditions, which results in the prediction 
of a steeper mass distribution in the power law range of 
masses in the HD case than in the MHD one, but no 
direct reference to instabilities in the layers. 

Figure 3 shows the compensated power spectra of the 
four main simulations. The power spectra are defined 
as the squared of the modulus of the Fourier transform 
of the velocity, integrated over a wave-number shell. If 
u.;(k) is the Fourier transform of the i velocity com- 
ponent, Uj(r), the power spectrum is Pi{k) = X)u«u*, 
where the sum is over all i and all wave-numbers k in 
the shell k < |k| < k + dk. P(k) is proportional to the 
contribution to the mean square velocity from all wave- 
numbers in the shell k < |k| < k + dk. 

The plots in Figure 3 have been arbitrarily shifted in 
the vertical direction. Deviations of more than a factor 
of two from a power law fit are found only at wavenum- 
bers k > 100, so the turbulence is roughly scale-free for 
almost two orders of magnitude in wavenumbers. The 
power law slopes depend somewhat on the exact range 
of wavenumbers used in the fit. We have chosen to mea- 
sure the power spectrum slope in the range 3 < k < 20, 
because larger wavenumbers are affected by the bottle- 
neck effect (e.g. Falkovich 1994; Dobler et al. 2003; Hau- 
gen & Brandenburg 2004). If the least square fits were 
extended up to k = 100, to estimate an effective power 
spectrum slope relevant for the turbulent fragmentation 
process, the slopes would be only slightly different. We 
get (5 — 1.9 and 2.0 from the Stagger code in the MHD 
and HD regimes respectively. The Enzo code in the HD 
regime gives (3 — 1.9, and Zeus in the MHD regime 
= 2.2. The corresponding values of the exponent of 
the power law range of the mass distribution of unstable 
cores are, according to the model of turbulent fragmen- 
tation, x = 1.4 and 3 for the MHD and HD regimes of 
the Stagger code, and x = 1.7 and 2.5 for the MHD and 
HD regimes of Zeus and Enzo respectively. 

The power spectra show that the Stagger code is only 
slightly more dissipative than Enzo, while Zeus is sig- 
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Fig. 7. — Test of convergence of the mass distribution with 
decreasing value of the density resolution parameter, /, from the 
Stagger-Code MHD run. The mass distribution is well converged 
at / = 8%. 



nificantly more dissipative than both. The numerical 
dissipation results in a power spectrum that is well fit 
by an extended power law, partly because an effect of 
the numerical dissipation is to suppress the bump of ex- 
cess power corresponding to the bottleneck effect, visi- 
ble especially in the least dissipative code, Enzo, around 
k = 40. The power law form of the power spectrum may 
therefore be deceiving (as evidence of numerical conver- 
gence), because it's slope is apparently dependent on the 
numerical diffusivity. It is only by comparing different 
codes or the same code with different values of the dif- 
fusivities that we can find a truly converged power spec- 
trum. If the power spectrum is not converged to its cor- 
rect slope, due to a lack of dynamic range of scales or 
an excess of numerical diffusivity, the slope of the mass 
distribution may be strongly affected. For example, in 
the MHD case we get x — 1.4 from the Stagger code, and 
x = 1.7 from the Zeus code. It is possible that even the 
Stagger code power spectrum is not fully converged, and 
that the correct value is approximately (3 — 1.8, in which 
case the slope is x = 1.36. This is suggested by the fact 
that in the HD regime PPM has a slightly smaller expo- 
nent than in the HD regime of the Stagger code. Padoan 
et al. (2006) have recently obtained an accurate mea- 
surement of the velocity power spectrum in the Perseus 
molecular cloud complex. Their result is [3 = 1.81 ±0.10, 
accurate enough to rule out the significantly larger power 
spectrum slopes produced by more dissipative SPH sim- 
ulations (see § 4). 

3. MASS DISTRIBUTIONS 

We compute the mass distribution of gravitationally 
unstable cores formed in turbulence simulations with- 
out self-gravity primarily to learn about the effect of 
turbulence under different conditions, for example with 
and without magnetic fields, and to compare with the 
predictions of the turbulent fragmentation model. The 
question of how self-gravity would modify the mass dis- 
tribution is a separate one, and will not be addressed 
in this work. However, it is important to stress that 
the present result are obtained after the driven turbu- 
lence has been statistically relaxed, which could not be 




0.1 1.0 10.0 

TO [M ] 



Fig. 8. — Same as Figure 7, but from the Zeus MHD run. 



achieved with self-gravity. The mass distributions de- 
rived in this work and the mass distribution predicted 
by Padoan & Nordlund (2002), can therefore represent, 
at best, a guess of the final outcome of more realistic sim- 
ulations with self-gravity. In such simulations including 
self-gravity the mass distribution of unstable cores may 
initially vary with time, as the most massive cores are 
still being assembled by converging turbulent flows while 
their central part has already collapsed. 

Cores are defined as connected overdensities that can- 
not be split into two or more overdensities of amplitude 
5p/ p > /. The unstable cores are simply the cores with 
mass larger than their Bonnor-Ebert mass. These defi- 
nitions are implemented in our clumpfind algorithm by 
scanning the density field with discrete density levels, 
each of amplitude / relative to the previous one. Only 
the connected regions above each density level that are 
larger than their Bonnor-Ebert mass are selected as un- 
stable cores. After this selection, the unstable cores from 
all levels form a hierarchy tree. Only the final (unsplit) 
core of each branch is retained. It is important to im- 
pose the Bonnor-Ebert condition while building the tree, 
otherwise some large unstable cores would be incorrectly 
eliminated if they split into smaller cores that were later 
rejected based on the Bonnor-Ebert condition. 

Clumpfind algorithms differ in the way they assign the 
surrounding mass to the cores. The popular algorithm by 
Williams, de Geus, & Blitz (1994), for example, uses up 
all the available mass (see their Figure 2). This results in 
a core formation efficiency of 100% above the threshold 
density, which is an artifact of that algorithm with no 
physical justification, though it may mimic a process of 
competitive accretion. Our algorithm, instead, assigns 
to each core only the mass within the density isosurface 
that defines the core (below that density level the core 
would be merged with its next neighbor). We prefer this 
definition because the smallest possible mass is assigned 
to each core and we want to study the effect of turbulent 
fragmentation in isolation, not trying to guess the out- 
come of the subsequent accretion. It turns out that, with 
this definition, and under conditions typical of molecular 
clouds, the unstable cores contain a few percent of the 
total mass. This suggests that not much of the remain- 
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Fig. 9. — Mean density distributions of the unstable cores se- 
lected from two snapshots of the Stagger Code MHD experiment 
(solid histogram), assuming a mean gas density of 10 4 cm -3 and 
a 6 pc size, as in previous figures. The dotted curve is the pdf 
of gas density of the whole computational domain, computed for 
the same two snapshots. The solid straight line shows the power 
law distribution of core mean densities that would be derived from 
the velocity scaling, N(n) oc n~ e /(P~ 4, using the specific value of 
P = 1.94 found in those two snapshots. The vertical dashed line 
marks the mean gas density. The two curves and the power law 
are normalized to approximately match each other at the largest 
densities, to illustrate that a very large fraction of the densest gas 
is found within unstable cores, while gas of decreasing density is 
increasingly harder to capture in unstable regions. 



ing mass will ever be accreted, as we know that the star 
formation efficiency in molecular clouds is approximately 
a few percent. 

Once the physical size and mean density of the system 
are chosen, the clumpfind algorithm depends only on two 
parameters: i) The spacing of the discrete density levels, 
/, and ii) the minimum density above which cores are 
selected, p m i n - In principle there is no need to define a 
minimum density, but in practice it speeds up the algo- 
rithm. We have verified that results (including the total 
mass in cores) do not change significantly for values of 
Pmin below the mean gas density, so we scan the density 
field only above the mean density. Notice that only half 
of the volume, but most of the mass, is found above the 
mean density, because, according to the lognormal pdf, 
most of the mass is packed in a small volume fraction. 

The parameter / may be chosen according to a phys- 
ical model providing the value of the smallest density 
fluctuation that could collapse separately from its con- 
tracting background. However, given the difficulty of 
predicting the outcome of the gravitational fragmenta- 
tion, we prefer to simply search for a convergence of the 
mass distribution with decreasing values of /. Luckily, 
the convergence is typically obtained already at a value 
of / w 16%, meaning that differences between the mass 
distributions with / = 8% and 16% are generally insignif- 
icant. The cores are therefore well defined, and in most 
cases they correspond to density fluctuations, relative to 
the surrounding gas, even much larger than /. 

In Figure 4 the mass distributions above 1 Mq are 
plotted for the main four experiments scaled to a mean 
density of 10 4 cm~ 3 , a box size of 6 pc, and a clumpfind 
density resolution / = 8%. Overplotted on the corre- 
sponding power law section of each mass distribution, 
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Fig. 10. — Same as Figure[9] but for one snapshot of the Stagger 
Code HD run. The power law here is N(n) oc n~ 3 /(' 3— 1 ) ! with the 
value of j3 = 2.06 computed for this snapshot. 



the dashed lines show the power law derived from the 
power spectrum slope and the shock jump conditions of 
each simulation, according to the turbulent fragmenta- 
tion model, x — 3/(4-/3) in the MHD regime, and 
x = 3/(5 — 2/3) in the HD regime. The general trend 
is recovered well, despite deviations to be expected be- 
cause this mass distributions are from single snapshots, 
not time averages. 

Figure 5 shows the mass distributions of the HD and 
MHD regimes (Enzo and Stagger Code respectively), 
computed with / = 16% and assuming a mean gas den- 
sity of 10 4 cm~ 3 . Each mass distribution is the result of 
matching two mass distributions, computed for a compu- 
tational box size of 1 pc and 6 pc. The 6 pc case makes 
it possible to sample masses in the range 1 — 10 Mq, 
and hence to probe the effect of the turbulence power 
spectrum and shock jump conditions on the mass distri- 
bution, but suffers from incompleteness for stars below 
1 Mq. The 1 pc case samples well the turnover region, 
and hence defines the peak mass for that mean density 
and rms Mach number, but does not yield intermediate 
and high mass stars. Figure 6 is equivalent to Figure 5, 
but uses the MHD Zeus run, instead of the Stagger Code 
run. 

The numerical mass distributions reproduce the sharp 
difference between the HD and MHD regimes predicted 
by the turbulent fragmentation model. The steeper mass 
distribution expected from the Zeus run, compared with 
the Stagger Code run, due to the steeper Zeus turbulence 
power spectrum, is also recovered. The slopes predicted 
by the turbulent fragmentation model are overplotted in 
each figure. Furthermore, the MHD regime yields a mass 
distribution of gravitationally unstable cores practically 
indistinguishable from Chabrier's stellar IMF (Chabrier 
2003), both in the Zeus and in the Stagger Code runs. 
Finally, the Stagger Code HD run also yields a mass dis- 
tribution consistent with the model prediction (see for 
example Figure 4). The relation between the mass dis- 
tribution and the power spectrum and shock jump condi- 
tions is therefore successfully tested with 3 different codes 
at very high numerical resolution. As discussed below, 
the model predictions for the HD regime are also qual- 
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Fig. 11. — Mean density versus mass, of the unstable cores se- 
lected from two snapshots of the Stagger Code MHD run, assuming 
a mean gas density of 10 4 cm -3 and a 6 pc size, as in previous fig- 
ures. The solid line on the left marks the Bonnor-Ebert mass, 
while the solid line on the right a mass 100 times larger than the 
Bonnor-Ebert value. The turbulence produces a range of masses 
and mean densities of unstable cores spanning three orders of mag- 
nitude. The dashed line is the power law relation between mean 
density and mass based on the average scale dependence of these 
two quantities, n oc W/C 8- 2 W ; using (3 = 1.94, as in Figure|9] 



itatively confirmed with a fourth code (the TVD code), 
at lower resolution. 

To verify the convergence of the clumpfind algorithm 
with respect to the density resolution, expressed by the 
parameter /, we plot in Figures 7 and 8 the mass distri- 
butions of the MHD regime assuming a mean gas density 
of 10 4 cm -3 and a 6 pc size. Figure 7 is from the Stag- 
ger Code MHD run and Figure 8 from the Zeus MHD 
run. For this convergence study we have computed to- 
gether the mass distributions of unstable cores selected 
from three different snapshots. With these larger sam- 
ples, statistical deviations are reduced, resulting in a 
more sensitive test of convergence. Between / = 32% 
and / = 8% there is a tendency to fragment the largest 
cores and create a larger number of small cores. However, 
the differences between / = 8% and / = 2% are rather 
small. Furthemore, the slope of the mass distribution 
above 2-3 Mq is rather independent of resolution, even 
if the power law section of the mass distribution tend to 
move slightly to larger masses at very large values of /. 

4. DISCUSSION 

4.1. Variance Versus Scaling 

The turbulent fragmentation model estimates the mass 
distribution of unstable cores, based on the assumption 
that their size is determined by the thickness of post- 
shock layers. So far we have used the simulations and the 
clumpfind algorithm to confirm the predicted relation be- 
tween the slope of the mass distribution and that of the 
velocity power spectrum. The model also implies that, 
on the average, the size and mean density of cores should 
be scale dependent as well, at least for cores significantly 
more massive than the average Bonnor-Ebert mass. In 
the absence of variance, the distributions of core proper- 
ties are derived from their scale dependence assuming the 
self-similar distribution of scales, N(L) oc L~ 3 . A quan- 
tity X, scaling like X oc L a , would have a power law 
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Fig. 12. — Same as Figure [TT] but for the cores selected from 
one snapshot of the Stagger Code HD run. The dashed power law 
here is n oc m 

(13-1) / (5-2/3) ^ with the vame of p _ 2.06 computed 
for this snapshot. 



distribution N(X) oc X~ 3 /°, where the distributions are 
expressed as probability densities per logarithmic inter- 
vals (or, equivalcntly, they can be thought of as cumula- 
tive functions derived from the integration of probability 
densities per linear intervals). A very steep distribution 
of a quantity X is therefore equivalent to a weak scale 
dependence of X (small values of the exponent a), re- 
sulting in a limited range of values of X caused by its 
scale dependence alone. 

However, the properties of cores arising in turbulent su- 
personic flows are random variables following some distri- 
butions that can be roughly characterized by their stan- 
dard deviation. This must be true on any scale. If the 
standard deviation of the distribution at a fixed scale is 
very large, it may exceed the range of values expected 
from the scale dependence, and the full distribution (in- 
tegrated over all scales) may resemble more the distribu- 
tion at a fixed scale than the power law predicted from 
the scale dependence. In the case of the core mass dis- 
tributions in both the MHD and HD regimes, a > 1, 
meaning that the logarithmic range of core masses ex- 
ceeds the logarithmic range of scales, and therefore the 
scale dependence is expected to leave a strong imprint 
in the mass distributions. This is a case in which the 
scale dependence of a quantity results in a well defined 
power law distribution of that quantity. If instead the 
scale dependence of a quantity X is weak (a < 1), it may 
be completely masked by the variance of its distribution 
at a fixed scale. 

Let's consider first the distribution of the core mean 
density. According to our model, the scaling of veloc- 
ity results in the core mean density distribution N(n) oc 
n -6/{0-i) _ n -6.67 ; in tne MHD re gi me; an d jy(n) ^ 

n -3/(/3-i) _ n -3.33 m |- ne jjQ re gi me; where the numeri- 
cal estimate of the exponents assumes f3 = 1.9. These 
density distributions are both steeper than the corre- 
sponding mass distributions, while the effect of the vari- 
ance in density should be comparable to that in mass, 
as mass is proportional to density. Therefore, even if the 
mass distributions have power law tails, the core mean 
density distributions may not show any power law. This 
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Fig. 13. — Size distribution of the unstable cores selected from 
the same two snapshots of the Stagger Code MHD experiment as 
in the left panels of Figures [9] and 1111 The core size is defined as 
the cubic root of its volume. The power law distribution resulting 
from the scale dependence of the core size, IS 
plotted as a solid line, using the value of f) = 1.94 found in the two 
snapshots. 



is certainly the case in the MHD regime, given the very 
large exponent and hence the very weak scale dependence 
of the mean density (very narrow range of densities over 
a large range of scales). 

In Figures [9] and [10] we show the MHD and HD mean 
density distributions of unstable cores from the Stag- 
ger Code experiments (solid histograms) , the power laws 
based on the velocity scaling (straight solid lines) and the 
pdf of gas density from the whole computational volume 
(dotted curves). Consistent with the arguments above, 
in the MHD regime there is no evidence of the power law 
scaling, while in the HD regime a short power law is vis- 
ible, with a slope consistent with the model prediction. 
Notice how the core mean density distribution gradually 
departs from the general density pdf toward lower den- 
sities, as a result of the selection of only gravitationally 
unstable cores. This large amount of gas mass not locked 
in unstable cores, even at relatively high densities, illus- 
trates how the turbulence controls the efficiency of the 
star formation process. 

The weak scale dependence of the core mean density 
can also be appreciated by plotting the core mean density 
versus the core mass, as done in Figures ITT1 and [T2"l These 
scatter plots don't show much evidence of the power law 
relation of mean density and mass from the scale depen- 
dence of these two quantities, n oc m^^ 1 ^^ 2 ^ = m 021 
in MHD, and n oc m (/3-i)/(5-2/3) = m o.75 in HDj as _ 

suming again (3 — 1.9 in the numerical estimates. At 
any given mass, the density may span almost three or- 
ders of magnitude, although most values are within an 
factor of ten range. This large variance is unavoidable 
in isothermal supersonic turbulence, due to the broad 
Lognormal density pdf, with standard deviation approx- 
imately equal to half the rms Mach number of the tur- 
bulence. 

The average value of the core size should also be 
scale dependent, with the corresponding distributions 
being N(l) cx l~ 6 /^-f 3 ) = r 5 45 in MHD and N(l) cx 
1-3/(2-/3) _ ;-30 m HD ( assum i ng p = 1.9 m tne nu _ 
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Fig. 14. — Same as in Figure [T3j but for the Stagger Code HD 
experiment. The scaling law would now give an extremely steep 
power law, N(l) oc i -3 /( 2- ' 3 ) ) or essentially no size range. 



merical estimates). The MHD power law is rather steep, 
but the density variance translates into a smaller size 
variance, as I oc n -1 / 3 , so even a rather weak scale de- 
pendence may result in a power law size distribution. In 
the HD regime, there is almost a perfect cancellation of 
the scale dependence of the shocked layer column den- 
sity and volume density, making the average thickness 
of such a layer almost scale independent. The resulting 
power law distribution would be so steep, that even a 
small variance would completely mask it. 

Figures [13] and [14] show the core size distribution of 
the same cores as in Figures [S] to fTS] The core size is 
defined as the cubic root of its volume. The power law 
distribution of the mean core size resulting from the scale 
dependence is plotted as a solid line in the MHD case, 
with the exponent corresponding to a value of (3 = 1.94, 
as measured from the two snapshots used to select the 
cores. In the MHD experiment there is some evidence of 
a power law tail, with the same slope as inferred from 
the scale dependence of our model. In the HD regime, 
instead, there is no evidence of the power law tail, as 
that would correspond to a sharp cutoff, which is easily 
masked by the variance of the size distribution at a fixed 
scale. 

In summary, based on the competition between vari- 
ance and scaling in establishing statistical distributions 
of core properties, we expect to find power law tails in the 
distributions of masses and sizes in the MHD regime, and 
only masses (or marginally densities) in the HD regime. 
The unstable cores selected from our numerical experi- 
ments have clear power law tails only in the case where 
these are expected, and in such cases the power law 
slopes are those predicted by the model. In Padoan and 
Nordlund (2002) the competition between the variance 
and the scaling is responsible for the flattening and the 
turnover of the mass distribution toward smaller masses, 
because toward smaller masses the selection of the unsta- 
ble cores becomes increasingly sensitive to their density. 
At large masses, because the scaling is the dominant ef- 
fect, the variance is neglected. 

In numerical simulations, even when the variance 
should not be dominant, power law tails may be absent 



9 




1 10 100 

k 



Fig. 15. — Power spectra compensated for the slope of the 
Stagger-Code HD run, /3 = 1.9. The TVD and SPH power spectra 
are the same as in Figure 2 of Ballesteros-Paredes et al. (2006), 
for the Mach numbers 3 and 6. 



from statistical distributions of core properties, as a re- 
sult of the limited range of scales relative to the actual 
interstellar turbulence. If the range of scales is reduced, 
the range of values in core properties due to the scale 
dependence is also reduced, possibly to the point of be- 
coming smaller than the variance of the distribution at a 
fixed scale. Lognormal-like tails may then cut short the 
power law distributions, as a numerical effect. 

It is important to appreciate that a driven turbulent 
flow may experience, over time, significant deviations 
from its average scaling laws, and that this may be the ex- 
planation for observed variations of the stellar IMF from 
place to place much in excess of the Poisson variance 
related to the statistical sample size. The scaling laws 
were understood phenomenologically by Kolmogorov as 
due to a scale independent energy dissipation rate, aris- 
ing from an efficient energy cascade from large to small 
scales in turbulent flows. This transfer from large to 
small scales takes approximately a dynamical time of the 
outer scale. Therefore, in a driven flow, any variations 
of the energy injection rate on a time-scale of order the 
dynamical time causes a "bump" in inertial-range scaling 
laws that has to propagate down the turbulent cascade, 
until it reaches the small viscous scales after approxi- 
mately a dynamical time of the outer scale. Because 
the typical lifetime of star-forming regions is compara- 
ble to this dynamical time (and star formation starts 
immediately when a molecular cloud is assembled), the 
turbulence can hardly be considered relaxed, and large 
variations of the IMF from place to place should be ex- 
pected. These variations should not be interpreted as the 
lack of a universal process of star formation, but rather 
as the evidence of both its turbulent origin and its short 
lifetime. 

4.2. Previous Results 

Ballesteros-Paredes et al. (2006) argue that the frag- 
mentation model of Padoan and Nordlund (2002) is in 
contradiction with their numerical results, based on TVD 
and SPH simulations without magnetic fields. They con- 
clude that the core mass distribution depends on the rms 
Mach number, but fail to point out that the Padoan and 
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Fig. 16. — Mass distributions of gravitationally unstable cores 
from the Enzo HD runs with M s = 6 and M s = 10, / = 8% 
and assuming a mean gas density of 10 4 cm -3 and a box size of 
6 pc. Each mass distribution contains unstable cores from two 
snapshots. The Chabrier (2003) IMF (dashed line) and the power 
law predicted by the turbulent fragmentation model (solid line) are 
also plotted. 



Nordlund model contains such a Mach number depen- 
dence, with the peak of the mass distribution shifting 
to lower masses as the Mach number increases, in agree- 
ment with the numerical results in Ballesteros-Paredes et 
al. (2006). In the Padoan and Nordlund model, the slope 
of the mass distribution for masses above the peak is in- 
dependent of the Mach number, also in agreement with 
the results of Ballesteros-Paredes et al. (2002) based on 
the TVD simulations (see their Figure 4), but in contra- 
diction with their SPH simulations (see their Figure 5). 

Figure [15] compares the power spectrum from the 
Stagger-Code HD run with two TVD and two SPH power 
spectra from Ballesteros-Paredes et al. (2006), for Mach 
numbers 3 and 6. The inertial range in both the TVD 
and SPH cases is not very extended, due to the low 
numerical resolution. The TVD code gives a slope of 
(3 ~ 2.2, the same value found in the Zeus run, for 
both Mach numbers. The extent of the inertial range 
in the TVD run is also comparable to the Zeus result 
at the same resolution (not shown). The power spec- 
tra of the SPH runs are instead much steeper, and their 
slope increases with decreasing Mach number, (3 w 2.7 for 
Ms fa 6 and » 2.9 for M s pa 3. As shown by the TVD 
runs, the power spectrum should not vary much with 
Mach number between Ms = 6 and Ms = 3. For lower 
Mach numbers, the power spectrum should become shal- 
lower, and converge to a value of (3 ~ 5/3 for Ms < 1. 
The SPH power spectrum slope is therefore much too 
steep and its Mach number dependence unphysical. 

In summary, it appears that the TVD runs of 
Ballesteros-Paredes et al. (2006) are able to qualitatively 
reproduce the turbulent fragmentation process that we 
have tested in the present work with much larger numer- 
ical resolution, with three different grid-based codes and 
both with and without magnetic fields. The complete ab- 
sence of an inertial range with a reasonable slope, or with 
a reasonable dependence of the slope on the Mach num- 
ber, makes their SPH simulations totally inadequate for 
testing the turbulent fragmentation model, as the model 
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relies on the scale-free nature of turbulent flows. Nev- 
ertheless, Ballesteros-Paredes et al. (2006) seem to base 
their conclusions partly on the SPH results, even if in 
contradiction with their own more robust TVD results 
that are confirmed here. 

Ballesteros-Paredes et al. (2006) try to argue that the 
magnetic field plays no role in the IMF, to justify the rele- 
vance of their simulations without magnetic fields. How- 
ever, we have shown here that the HD regime produces 
much steeper mass distributions than the MHD regime. 
The HD regime is therefore a poor choice of physical pa- 
rameters, if the aim is to extract the power law mass 
distribution from simulations with a very limited numer- 
ical resolution, as it is much more variance-dominated 
than the MHD regime. As explained above, it is not 
surprising that the very limited tails of their power laws 
may possibly appear more Lognormal than straight. 

To address directly the issue of the Mach number de- 
pendence of the mass distribution, raised by Balesteros- 
Paredes et al. (2006), we plot in Figure [T()] the mass 
distributions from the Enzo HD runs with M§ = 6 and 
Ms = 10 and with / = 8%. These mass distributions 
are computed assuming a mean gas density of 10 4 cm~ 3 
and a box size of 6 pc. Each distribution contains cores 
from two snapshots. The Figure shows that the power 
law part of the mass distribution, above 1-2 Mq , is inde- 
pendent of the Mach number and matches the prediction 
of the turbulent fragmentation model, that is fc~ for 
the power spectrum slope f3 = 1.9 of the HD Enzo runs. 

5. CONCLUSIONS 

We have used large numerical simulations of supersonic 
MHD and HD turbulence to test the turbulent fragmen- 
tation model of Padoan and Nordlund (2002). The model 
predicts a power law distribution for large masses, related 
to the turbulence energy power spectrum slope, and the 
shock jump conditions. This power law mass distribu- 
tion is confirmed by the numerical experiments. The 
model also predicts that the HD regime should yield a 
much steeper mass distribution of unstable cores than 
the MHD regime, which is confirmed by the simulations. 
This feature of the fragmentation model is very impor- 
tant because it shows that even rather weak magnetic 
fields (super-Alfvenic turbulence) can be crucial in set- 
ting the initial conditions for the process of star forma- 
tion and in shaping the stellar IMF. 

While present-day star formation takes place probably 



always in the MHD regime, star formation at very high 
redshift may well occur in the HD regime, both because 
the field strength is still low and because the value of 
the critical magnetic field strength that defines the HD 
regime is larger at higher temperatures. The IMF of 
the earliest population II stars, and perhaps the latest 
population III stars as well, may be formed in turbulent 
environments in this HD regime, resulting in an IMF 
narrowly peaked around a mass of order 10 M©. The 
effect of such a peculiar IMF of early stellar populations 
on the ionization history and metallicity evolution of the 
universe should be investigated. 

Numerical simulations can quantitatively account for 
the fundamental role of the turbulence in setting the ini- 
tial conditions for the process of star formation only if 
they can generate an inertial range of turbulence, which 
requires both low numerical diffusivity and large numer- 
ical resolution. Furthermore, to model present-day star 
formation that occurs in the MHD regime, the magnetic 
field cannot be neglected, even if the turbulence is as- 
sumed to be super-Alfvenic. SPH simulations of large 
scale star formation to date fail in all three fronts: nu- 
merical diffusivity, numerical resolution and presence of 
magnetic fields. This should cast serious doubts on the 
value of comparing predictions based on SPH simulations 
with observational data (see also Agertz et al. 2006). 

Finally, the mass distribution of unstable cores found 
in the MHD simulations is indistinguishable from the 
Chabrier stellar IMF (Chabrier 2003) and in qualitative 
agreement with the less well determined mass distribu- 
tions of prestellar cores selected from dust emission or 
molecular line observations. Such a coincidence may in- 
dicate that gravitational fragmentation, competitive ac- 
cretion or merging, all absent in these turbulence simula- 
tions, may not play a major role in the origin of the stellar 
IMF, a fascinating idea to be tested with the next gen- 
eration of high resolution simulations of self-gravitating 
turbulence. 
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